Classical simulation of infinite-size quantum lattice systems in one spatial dimension 
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Invariance under translation is exploited to efficiently simulate one-dimensional quantum lattice 
systems in the limit of an infinite lattice. Both the computation of the ground state and the 
simulation of time evolution are considered. 

PACS numbers: 



Numerical renormalization group (RG) methods ob- 
tained a first remarkable success with Wilson's solution 
of the Kondo problem [l[ and, ever since the irruption of 
White's density matrix renormalization group (DMRG) 
algorithm Q, are solidly established as the dominant 
computational approach to quantum lattice systems in 
one spatial dimension (ID) |3j. Recently, progress in our 
understanding of quantum entanglement has prompted a 
series of new developments [1, [B|, H, 0, S, that extend 
the domain of RG methods in two main directions. On 
the one hand, still in the context of ID lattice systems, 
the time-evolving block decimation algorithm (TEBD) al- 
lows for the simulation of time evolution [4j. On the 
other, Verstraete and Cirac [5| have introduced projected 
entangled-pair states (PEPS) to address the simulation 
of quantum lattice systems in two and higher spatial di- 
mensions. 

In this work we propose an algorithm, based on TEBD, 
to simulate ID quantum lattice systems in the thermody- 
namic limit. Bulk properties of matter are best studied 
in an infinite system, where they are not contaminated 
by finite-size corrections or boundary effects. However, 
for most algorithms the cost of a simulation grows with 
the system size, and the thermodynamic limit can only 
be reached by extrapolating results for increasingly large 
systems. Here we exploit two facts, namely invariance 
under translations of the system and parallelizability of 
local updates in TEBD, to obtain a noticeably simple and 
fast algorithm, referred to as infinite TEBD or iTEBD, 
to simulate infinite systems directly, without resorting 
to costly, less accurate extrapolations. We describe the 
iTEBD algorithm and test its performance by computing 
the ground state and time evolution for a quantum spin 
chain in the thermodynamic limit. In addition to offer- 
ing a very competitive alternative to DMRG for infinite 
systems, iTEBD plays a key role in entanglement renor- 
malization techniques [lOfland in the extension of PEPS 
d to infinite 2D lattices 

We consider an infinite array of sites in ID, where each 
site r, r 6 Z, is described by a complex vector space 
y[ r ] ^ f g n i^ c dimension d. Let vector 



denote a pure state of the lattice and operator H, 

H = ^h [r ' r+1 \ (2) 

r 

a Hamiltonian with nearest neighbor interactions, and let 
and H be invariant under shifts by one lattice site 
12J. Given an initial state \^o), our goal is to simulate 
an evolution according to H , both in real time 



\* t ) = exp{-iHt)\* ), 
and in imaginary time (l3| . 

exp(-i2Y)|* ) 
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The TEBD algorithm represents |\&) through a matrix 
product state (MPS) [13]. See 0] for details on this struc- 
ture, that we briefly review for an infinite ID lattice. Let 
[<r] and [r+l>] denote the semi-infinite sublattices made 
of sites {— 00, • • • , r} and {r +1, • • • , 00}. If the Schmidt 
decomposition of |^) according to this bipartition reads 



(5) 



where we assume the Schmidt rank \ to be finite, then 
the spectral decomposition of the reduced density matri- 
ces for [<r] and [r +lD>] are 



,[<r] 



= E (A") 2 i<fL <rl x$L <rl 



(6) 



a=l 

p [r+i» = J- (AM) 2 |$L r+1>1 )($L r+1>1 |. (7) 

We use a three-index tensor rM to relate the Schmidt 
bases for two left (respectively right) sublattices: 

X d 



(i) 



= j;EWi*m (8) 

/3=1 <=1 

i$L r>I ) = EE r S lA r ] K w >K r+1>1 >- ( 9 ) 

/3=1 <=1 

In particular, can be expanded in the local basis 
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FIG. 1: Diagrammatic representation of several decomposi- 
tions of |\&) in terms of a network of tensors, see for more 
details. Each filled circle represents a tensor F and each edge 
represents an index, such as a or i. Hanging from an open 
side of an edge there is an (omitted) orthonormal basis, such 



as |<3?!^ ') or An empty diamond on top of an edge 

represents a set of weights A (Schmidt coefficients) for the 
corresponding index, (i) Schmidt decomposition, Eq. (T5J), ac- 
cording to semi-infinite sublattices [<0] and [1>]. The bases 
|$L <0] > and |$L 1>] > are linked through index a with weights 
\a\ (ii) Expansion of [St) in terms of the Schmidt bases for 
the semi-infinite sublattices [<0] and [2>] and an orthonor- 
mal basis for site 1, see Eq. (|10[> . (Hi) Expansion of in 
terms of orthonormal sets of vectors for sublattices [<0] and 
[3>] and sites 1 and 2. 



:M) for site r and in terms of A^r^A^ 1 ] as 



< r ]\ 



l*L r+1>1 ), (io) 



a,f}=l 1=1 



or for sites {r, r+1} in terms of X^T^ X^r^ , 
and so on, see Fig- CU) ■ We also recall that in the TEBD 
algorithm the evolution operator exp(— iHt) in Eq. (151) 
is expanded through a Suzuki- Trotter decomposition [151 ] 
as a sequence of small two-site gates 

jjlr.r+l] _ exp (_ ih [r,r+l] St ^ ft < 1, (11) 

which we arrange into gates U AB and U BA , 

U AB = (S?) U^ 2r+1 \ U BA = (9) U^- 1 ^. (12) 



Because state is shift invariant, it could be rep- 
resented with a MPS where and A^ are indepen- 
dent of r. However, we will partially break translational 
symmetry to simulate the action of gates (fl2"|) on \^>). 
Accordingly, we choose a MPS of the form 



p[2r] = r A } 
p[2r+l] _ pB 



X [2r] = X A^ 
X [2r+1] = X B 



r G Z. 



(13) 



As in the finite n case Q, the simulation of the time 
evolution, see Eq. ([3]), is achieved by updating the MPS 
so as to account for the repeated application of gates U AB 
and U BA on |^}. But for n — oo, the action of the gates 
preserve the invariance of the evolved state under shifts 
by two sites, see Fig. ©, and only tensors T A ,T B , X A 
and X B need to be updated - a task that is achieved 
through simple matrix manipulations, see Fig. ([3]). In 
other words, whereas for n sites the TEBD algorithm 
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FIG. 2: Two diagrammatical representations (see also Fig. 
(ft]) and ref. j8|]) for state U AB \fy): (i) tensor network repre- 
sentation containing both an MPS for and two-site gates 
U acting on each pair of sites {2r, 2r + 1}, Vr £ Z; (ii) new 
MPS for U AB \$!). Notice that both structures are invariant 
under shifts by two lattice sites and are completely specified 
by a small number of tensors, in spite of the fact that they 
represent a state of an infinite ID lattice. 
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FIG. 3: In order to update the MPS after gate U has been 
applied, see Fig. @, we first contract the tensor network (i) 
into a single tensor ® a ij~i (ii)- We then compute the singular 
value decomposition of according to the index bipartition 
[cm] : [j-y], namely 9 = Y^f} X \oti}0>4 Y l3\j-y\ as in ("()• We 
introduce \ B back into the network (iv) and form tensors 
f A and f in (iv) by attaching to X and Y the inverse of 
the Schmidt coefficients A s . All such matrix manipulations 
are achieved with 0(d 2 \ 2 ) space and 0(d 3 \ 3 ) and need to 
be performed only once in order to update the MPS for the 
whole infinite chain. 



requires 0(ndx 2 ) space to store an MPS and 0(nd 3 x 3 ) 
time to simulate a small evolution exp(-iHSt) for 
n = oo sites the iTEBD requires computational space and 
time that scale just as 0(d 2 x 2 ) and 0(d 3 \ 3 )- Key to such 
dramatic cost reduction by a factor n is the fact that, in 
contrast to other approaches [6j, we use a MPS based on 
the Schmidt decomposition, allowing for a parallelized, 
local update of tensors T and A. 

Finally, evolution in imaginary time, Eq. is also 
simulated with iTEBD by simply replacing the two-site 
unitary gates exp(— ih 5t) in Eq. (jTTJ) with non-unitary 
gates exp(— h St), St 1 [hH - 

We have tested the performance of iTEBD by com- 
puting the ground state and by simulating time evolution 
for the ID quantum Ising chain with transverse magnetic 
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FIG. 4: Left: spectrum p a = (A Q ) 2 of the density matrix 
for half of the infinite chain, see Eq. ([SJl, for different val- 
ues of the magnetic field h z . For h z 3> 1 the eigenvalues p a 
decay very fast with a. As h z approaches the critical value 
ft* = 1, the decay is less pronounced and the computation of 
the spectrum becomes more demanding. Right: however, for 
h z = 1.001 we still obtain the 50, 100 and 250 first eigenval- 
ues with remarkable accuracy, the error on p a being typically 
several orders of magnitude smaller that its value. 

field, as defined by the Hamiltonian 

H = J2(4 ] 4 +1] + h*4 ] )- (14) 

This model is exactly solvable [l7j], allowing for a direct 
comparison of numerical results with the exact solution. 

Firstly, by simulating an evolution in imaginary time, 
Eq. ([3]), we have obtained an approximate MPS repre- 
sentation of the ground state \^> g ) for several values of 
h z > 1. We assess the accuracy of the result by focussing 
on the following: (i) the spectrum p a — (Aa') 2 of the re- 
duced density matrix for half the infinite lattice, Eq. 
which characterizes entanglement accross the boundary 
of two halves of the chain; and (ii) the two-point corre- 
lator 

C 2 (r) = <* g |4 0] 4 rl l* ff > ~ «* 9 |4 01 |* S » 2 , (15) 

which quantifies the strength of correlations between two 
spins that are r lattice sites apart. 

In both cases, we obtain quantitative agreement with 
the exact solution up to several significant digits, see 
Figs. (01 and |J5]). Computing the ground state |\& g ) 
is particularly fast and accurate for large values of h z 
(say h z > 1.1) and, as expected, becomes more resource 
intensive as h z approaches the critical value h* = 1 [l||. 
We find, however, that for values very close to the crit- 
ical point, such as h z = 1.001, we still obtain accurate 
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FIG. 5: Two point correlator C 2 (r) of Eq. |15J). Up: for 
h z = 1.001, when the distance r between spins is up to several 
hundred sites, Ci{r) seems to decay as a power law (as in the 
critical case) but for r of the order of 1000 the exponential 
character of the decay becomes manifest. Down: even though 
iTEBD is unable to properly reproduce the spectrum p a of 
half an infinite chain at criticality, h z = 1, (see Fig. Q), we 
still obtain a very accurate approximation to the correlator 
Gi{f) when r is of the order of several thousands of spins. 

approximations to the half-chain spectrum {p a } and the 
correlator C^ir) through affordable simulations, see Fig. 
([S]). And, most remarkably, reliable results for C*2(r) are 
obtained even at the critical point h z = 1 by using a 
MPS with a reasonably small \i m spite of the fact that 
such a MPS is no longer able to reproduce {p a }- 

Secondly, by switching to real time, we have simu- 
lated the evolution of the infinite system, prepared in 
the ground state of H in Eq. (| 14[) for h z = 10, when the 
magnetic field is abruptly modified from its initial value 
to h' z = 3. Fig. (J6j> shows the evolution in time of the 
magnetization 

m z (t) = (* g \4 ] \* 9 ), (16) 

which after several oscillations becomes stable at some 
value different from the initial one. 

Summarizing, with modest computational resources 
the iTEBD is able to analyze an infinite ID system, 
near to and at a quantum critical point, for which it 
obtained accurate results for quantities involving up to 
several thousands of sites. Comparable results with pre- 
vious approaches 0, 0> S| would require computations 
that are, to the best of our knowledge, unaffordable. In- 
deed, simulating a ID lattice of, say, n = 10, 000 sites, 
would have a time and memory cost proportional to n, 
and would still include 0(l/n) errors due to the finite size 
of the system, in addition to more important deviations 
near the boundaries of the lattice. 
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[12] Simple extensions of the algorithm include: (i) interac- 
tions /i' r,s l with longer range, \r — s\ > 1, (ii) interactions 
involving more than just two sites (e.g. /i' r,3,t '"'), (Hi) 
time-dependent hamiltonians and (iv) systems invariant 
under shifts by m sites, with m > 1. 

[13] For a finite lattice with n sites and a hamiltonian H with 
gap A > 0, simulating imaginary time evolution as in 
Eq. (j4| for large r yields a good approximation \^f r ) to 
the ground state |vl/ 9 ) of H , since one can show that 



time [1/h] 



|(*r|* 9 )| >l-0(- 



l(*o|* 9 )|. 



(17) 



FIG. 6: Simulation of the evolution in time of m z (t) in Eq. 
(|16l) . The infinite chain is initially in the ground state of H 
with h z = 10. At t = the magnetic field is changed to 
h z — 3. Up: the magnetization, originally at a value m z (0) « 
0.5, oscillates in time until it stabilizes around a lower value 
m z (oo) ~ 0.48. By considering increasingly large values of x 
we are able to simulate the evolution of m z (t) for long times 
with remarkable accuracy. 



Work in progress considers a range of applications of 
the iTEBD, including the characterization of quantum 
critical points and cross-over between quantum phases, 
and the study of the response of ID systems to external 
probes. Originally developed to help in the context of 
entanglement renormalization , the method is also the 
key to extend 2D PEPS 
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to infinite systems 
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